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We study by extensive Monte Carlo simulations the transport of itinerant spins travelling inside 
a multilayer composed of three ferromagnetic films antiferromagnetically coupled to each other in a 
sandwich structure. The two exterior films interact with the middle one through non magnetic spac- 
ers. The spin model is the Ising one and the in-plane transport is considered. Various interactions 
are taken into account. We show that the current of the itinerant spins going through this system 
depends strongly on the magnetic ordering of the multilayer: at temperatures T below (above) the 
transition temperature Tc, a strong (weak) current is observed. This results in a strong jump of the 
resistance across Tc- Moreover, we observe an anomalous variation, namely a peak, of the spin cur- 
rent in the critical region just above Tc. We show that this peak is due to the formation of domains 
in the temperature region between the low-T ordered phase and the true paramagnetic disordered 
phase. The existence of such domains is known in the theory of critical phenomena. The behavior of 
the resistance obtained here is compared to a recent experiment. An excellent agreement with our 
physical interpretation is observed. We also show and discuss effects of various physical parameters 
entering our model such as interaction range, strength of electric and magnetic fields and magnetic 
film and non magnetic spacer thicknesses. 

PACS numbers: 72.25.-b, 73. 21. Ac, 75.75.-|-a 



I. INTRODUCTION 

The so-called giant magnctorcsistance (GMR) was dis- 
covered experimentally twenty years ago by Fert and 
coworkersi and by Grunberg^. Since then, intensive 
investigations, both experimentally and theoretically, 
have been carried out to understand the origin and the 
behaviors of the spin current in magnetic multilayer 
systems i^i^'^'i'^'^ This spectacular development is due 
mainly to an important number of industrial applica- 
tions using such systems in data storage and magnetic 



sensors 



3.4 



Experimental observations show that when the spin of 
an itinerant spin is parallel to the spins of the environ- 
ment it will go through easily while it will be stopped 
if it encounters an antiparallel spin medium. The resis- 
tance is stronger in the latter case resulting in a GMR. 
Although many theoretical investigations have been car- 
ried out, detailed understanding of the influence of each 
physical parameter on the spin current is still lacking. For 
example the role of interface scattering and the effect of 
interface roughness on the GMR are still currently inves- 
tigated (see for example RefsH0 and references therein). 
In addition, to date no Monte Carlo (MC) simulations 
have been performed regarding the temperature depen- 
dence of the dynamics of spins participating in the cur- 
rent. This defines the aim of this work. 

This paper deals with the transport of spins in a sys- 
tem composed of three magnetic films. We show that 
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the spin current depends on the orientation of the lattice 
spins found on the trajectory. The dependence of the spin 
transport on the magnetic ordering, i.e., on the tempera- 
ture is studied. The difhculty of the task is that we have 
to deal at the same time with surface and interface effects 
and with dynamical properties of itinerant spins interact- 
ing with the lattice spins. Surface physics of systems such 
as films and multilayers have been enormously studied 
at equilibrium during the last 30 years. This was moti- 
vated in particular by applications in magnetic recording, 
let alone fundamental theoretical interests. Much is un- 
derstood theoretically and experimentally in thin films 
whose surfaces are 'clean', i.e., contain no impurities, no 
steps etc i^^i-'^^'^^i-'^^'^'^i^^ Far less is known — at least theo- 
retically — on complicated thin films with special surface 
conditions such as defects, arrays of dots and magneti- 
zation reversal phenomena. As a result, studying the 
behavior of itinerant electrons injected into such systems 
is a formidable task which cannot be fulfilled in every 
respect. 

The paper is organized as follows. Section II is devoted 
to the description of our model and the rules that gov- 
ern its dynamics. We take into account (i) interactions 
between itinerant and lattice spins, and (ii) interactions 
between itinerant spins themselves (iii) interactions be- 
tween lattice spins. Where rules governing the dynamics 
are concerned, we include a thermodynamic force due to 
the gradient of itinerant spin concentration, an applied 
electric field that drives electrons, and the effect of a mag- 
netic field. In section III, we describe our MC method 
and discuss the results we obtained for several physical 
quantities in various situations, e.g., the mean free-path, 
the spin current and the resistance. Comparison with 
a very recent experimentii^ performed on a permalloy- 
insulator multilayer is also shown in this section. Con- 
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eluding remarks are given in Section IV. 



that the strength of Jjj- is physically reasonable, 
shape of the interaction is sketched in Fig. [2l 



The 



II. MODEL 



A. Interactions 



We consider in this paper three ferromagnetic films be- 
ing antiferromagnetically coupled to each other via non- 
magnetic layers. We use the Ising model and the face- 
centered cubic (FCC) lattice for the films. The system is 
shown in Fig. [1] where the films are stacked along the z 
direction. 

The multilayer is made up of three films each of which 
has a volume given by N^^ x Ny x N^-, where denotes 
the number of atomic layers (i.e., single film thickness). 
Periodic boundary conditions are used in the xy planes. 
Nonmagnetic spacers sandwiched between films have a 
thickness d. 



FIG. 1; Ground state spin configuration. Thick arrows indi- 
cate the spin orientations. Nonmagnetic spacers are black. 

Spins at FCC lattice sites are called "lattice spins" 
hereafter. They interact with each other through the 
following Hamiltonian: 



J' 



(1) 



where S,; is the Ising spin at lattice site i, indi- 
cates the sum over every nearest-neighbor (NN) spin pair 
(Si , Sj ) . For simplicity, we will consider the case where 
all exchange interactions j are ferromagnetic and equal 
to J(> 0), except for the interaction across the nonmag- 
netic spacer which we define using the following RKKY 
model, 



cos{arij) 



(2) 




FIG. 2: RKKY interaction for two spins across the non- 
magnetic layer is shown as a function of their distance rij. 
Jo = 16.6 and a — 6.18 have been used so that Jij=-1.005 
for nonmagnetic spacer thickness d = 2. 

When the coupling across nonmagnetic layers is anti- 
ferromagnetic, the ground state corresponds to the two 
exterior films having spins pointing in one direction and 
the interior one spins pointing in the opposite direction. 

In order to study the spin transport inside the multi- 
layer system described above, we consider a flow of itiner- 
ant spins interacting with each other and with the lattice 
spins. The interaction between itinerant spins is defined 
as follows, 



Tin 



(3) 



where is the Ising spin at position and j) 
notes a sum over every spin pair (s^, Sj). The interaction 
Ki j depends on the distance between the two spins, i.e.. 



n, = 



- rjl- A specific form of Ki j will be chosen be- 
low. The interaction between itinerant spins and lattice 
spins is given by 



(id) 



(4) 



where the interaction lij depends on the distance be- 
tween the itinerant spin and the lattice spin S^. For 
the sake of simplicity, we assume the same form for K^ j 
and li j , namely. 



= A'oexp(-ry) 
lij = /ocxp(-r,j). 



(5) 
(6) 



where Kq and Iq are constants expressing the respective 
strength of interactions. 



B. Dynamics 



Here, i and j refer to spins on either side of a nonmagnetic 
layer, and Jo a-nd a are constants chosen in such a way 



Let us now explain the procedure we utilize in our sim- 
ulation. First we study the thermodynamic properties of 
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the multilayer system alone, i.e., without itinerant spins, 
using Eq. ([1]). In this view, we perform MC simulations 
in order to determine quantities as the internal energy, 
the specific heat, layer magnetizations, the susceptibility, 
... as functions of temperature TJ^ From these physical 
quantities we determine the critical temperature Tc be- 
low which the system is in the ordered phase, e.g., with 
up-spin phase for the outer films and down-spin phase 
for the middle film. The total staggered lattice magne- 
tization is defined as M = {Mi - M2 + Ms)/^ where 
Mi{i = 1,2,3) is the magnetization of the ith film. We 
depict in Fig|3] the lattice magnetization versus T. 
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FIG. 3: Total staggered lattice magnetization versus temper- 
ature T. Mo is the ground state staggered lattice magnetiza- 
tion. Tc is ~ 9.75 in unit of J = 1. 

Figure |4] shows the susceptibility calculated from the 
fluctuations of A/ for two spacer thicknesses d ~ 2,3. 
Tc is equal to ~ 9.75 and ~ 9.49 respectively for these 
spacers. 




FIG. 4: Susceptibility x of staggered lattice magnetization 
versus temperature T for two spacer thicknesses d — 2 (solid 
line) and d = 3 (dashed line). 

Once the lattice has been equilibrated at T, we study 
the dynamics of itinerant spins at that temperature by in- 
jecting itinerant spins with a density n into the multilayer 
system. There are two ways of doing this: i) the itiner- 
ant spins move parallel to the film surface (GIF case); ii) 
the itinerant spins move perpendicular to the films (CPP 
case). In this paper we show results in the GIF case. 

The itinerant spins in the GIF case move into the sys- 
tem at one end, travel in the x direction, escape the sys- 
tem at the other end to reenter again at the first end 
under periodic boundary conditions (PBC). Note that 



PBC are used to ensure that the average density of itin- 
erant spins remains constant during the time (stationary 
regime). The dynamics of itinerant spins is governed by 
the following interactions: 

i) an electric field E is applied in the x direction. Its 
energy is given by 



He = -E- V, 



(7) 



where v is the velocity of the itinerant spin 

ii) a chemical potential term which depends on the 
concentration of itinerant spins within a sphere of radius 
D2 ("concentration gradient" effect). Its form is given 
by 



(8) 



where ?T.(r) is the concentration of itinerant spins in a 
sphere of radius D2 centered at r. is a constant taken 
equal to Kq for simplicity. 

iii) interactions between a given itinerant spin and lat- 
tice spins inside a sphere of radius Di (Eq. S]) . 

iv) interactions between a given itinerant spin and 
other itinerant spins inside a sphere of radius D2 (Eq.E]). 

Let us first consider the case without an applied mag- 
netic field. 

The simulation is carried out as follows: at a given T 
we calculate the energy of an itinerant spin by taking into 
account all the interactions described above. Then we 
tentatively move the spin under consideration to a new 
position with a step of length vq in an arbitrary direction. 
Note that this move is immediately rejected if the new 
position is inside a sphere of radius ro centered at a lattice 
spin or an itinerant spin. This excluded space emulates 
the Pauli exclusion principle in the one hand, and the 
interaction with lattice phonons on the other hand. For 
the sake of example, if the spacing between NN lattice 
spins is \/2 then ro is taken of the order of 0.05. This 
value can be made temperature-dependent to account for 
the cross section of phonon-electrons collisions. 

If the new position does not lie in a forbidden region of 
space, then the move is accepted with a probability given 
by the standard Metropolis algorithmic ; in particular, 
it is always accepted if the energy of the itinerant spin at 
the new position is lower than its previous value. 



III. MONTE CARLO RESULTS 

In this paragraph, we show the results obtained by MC 
simulations with the Hamiltonians given above. All Ising 
spins arc of magnitude s = S = I. 

The parameters we use in most calculations are, except 
otherwise stated, = 36, Ny = 10 and = 5 for the 
dimension of the films, d = 2 for the spacer thickness. 
We also make use of FBG in the xy plane. 

At each temperature the equilibration time for the lat- 
tice spins lies around 10^ MG steps per spin and we 
compute statistical averages over 10^ MC steps per spin. 
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FIG. 5: Mean free path A in unit of the FCC cell length versus FIG. 6: Mean free path A versus Di, at T = 1, D2 = 1 and 
temperature T , for 100 MC steps. E = 1. 



Taking J = 1. wc obtain Tc ~ 9.75 for the estimate of 
the critical temperature of the lattice spins (see FigslS] 
and 2]) . Before calculating the mean free path and the 
spin current, we let n itinerant spins travel through the 
system several thousands times until a steady state is 
reached. The parameters used for the results shown be- 
low are Di = D2 = 1 (in unit of the FCC cell length), 
Ko = /o = 2, n = 1500, vq ^ 1, tq = 0.05. 

In Fig. [5] we sketch the travelling length A computed 
after a fixed lapse of time as a function of temperature T. 
As can be seen, A is very large at T < Tc- We note that 
there is a small depression in the transition region. We 
will show below that this has important consequences on 
the spin current. We also note that at very low T (T < 4), 
the mean free path suffers a decrease with decreasing 
T. This is a well-known artefact of MC simulation at 
very low T: the moving probability is so small that the 
motion of itinerant spins is somewhat slowed down. As 
we will show below when comparing with experimental 
data, this freezing is also observed in real systems due to 
finite experimental observation time. 

Figures [S] and [7] show the effects of varying Di and 
D2 at a low temperature T = 1. As seen in Fig. [6l A 
is very large at small Di: this can be explained by the 
fact that for such small Di, itinerant spins do not "see" 
lattice spins in their interaction sphere so they move al- 
most in an empty space. The effect of D2 is on the other 
hand qualitatively very different from that of Di as seen 
in Fig. [7J A is saturated at small D2 and decreases to 
the minimum value, namely A=l, at large D2. We con- 
clude that both Di and D2 dominate A at their small 
values. However, at large values, only D2 has a strong 
effect on A. This effect comes naturally from the criterion 
on the itinerant spins concentration used in the moving 
procedure. 

The mean free path is shown in Fig. [8] as a function 
of T for two magnetic film thicknesses. In the absence of 
interface impurities, it is expected that there would be 
no large effects on the motion of itinerant spins. This 
is indeed what we observe here. Note however that the 
mean free path for the smaller magnetic film thickness is 
systematically smaller than that of the thicker film. We 
will discuss on the role of interfaces below while showing 
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FIG. 7: Mean free path A versus 1)2, at T = 1, Di = 1 and 
E = \. 

the resistance (Fig. [TT|) . 

We show in Fig. [9] the effect of the spacer thickness on 
the mean free path. Note that for each thickness value we 
have used the inter film coupling constant Ji calculated 
by Eq. [2] Increasing the thickness, i. e. decreasing J^, 
will result in a decrease of the mean free path visible at 
low T as can be seen in Fig. [51 This is expected since the 
itinerant spins at magnetic-nonmagnetic interfaces have 
weaker inter film coupling energy, so they are scattered 
more easily. 

We show in Fig. [TU]the effect of the electric field E for 
T both above and below T^- The low- field part verifies 
the Ohm regime. 
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FIG. 8: Mean free path A versus T, for several thickness values 
of the magnetic film, with D\ = D2 = 1 and E = 1. 
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FIG. 9: Mean free path A versus T, for several spacer thick- 
nesses with D\ — D2 — 1 and E — 1. 
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FIG. 11: Resistance R in arbitrary unit versus temperature 
T for two magnetic layer thicknesses. 





FIG. 10: Mean free path A versus E, below and above Tc 
with Di = D2 = 1. 



For the i-th layer, we define the resistivity as 

1 



(9) 



where is the number of spins crossing a unit area per- 
pendicular to the X direction per unit of time. Note that 
this definition is apphed to three magnetic {i = 1,3,5) 
and two nonmagnetic layers {i = 2,4). The total resis- 
tance R is defined as 



^ 1 



(10) 



This definition is suitable for low-T phase where the spin 
current is distinct in magnitude between magnetic and 
nonmagnetic layers. On the contrary, in the paramag- 
netic phase the spin current is almost spatially uniform, 



and the resistance can be defined as 

5 



(11) 



In Fig. [TT] we show the resistance i? as a function of 
temperature. 

There are several striking points: 

• i? is very low in the ordered phase and large in the 
paramagnetic phase. Below the transition temper- 
ature, there exists a single large cluster with small- 
sized excitation inside it (see Fig. [T3|) . so that any 
itinerant spin having the appropriate orientation 
goes through the structure without any hindrance. 
The resistance is thus very small. 

• R exhibits a cusp at the transition temperature, 
the existence of which was at first very surprising. 
While searching for its physical origin, we found 
that it was due to changes in the size distribution 
of clusters of lattice spins as the transition tem- 
perature is approached ; it is known indeed from 
the theory of critical phenomena that clusters of 
up (resp. down) spins of every size form when T 
approaches Tc in a critical phase transition. At 
Tc, the distribution of cluster sizes displays clus- 
ters of various sizes, as can be seen from Fig. [12] 
(more details on the cluster construction algorithm 
will be given below). As a result, the conductiv- 
ity is drastically lower than in the ordered phase 
since itinerant electrons have to steer around large 
clusters in order to go through the entire struc- 
ture. Indeed thermal fluctuations are still not large 
enough to allow the itinerant spin to overcome the 
energy barrier created by the opposite orientation 
of the clusters; this is all the more influential that 
wc flxed an infinite spin-fiip time, and this forbids 
the itinerant electron to reverse its orientation in 
order to reduce energy barriers. 

• Below there is no effect of magnetic layer thick- 
ness on R. However, for T > Tc, the larger thick- 
ness yields a smaller R. This can be explained 
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by the effect of interfaces at nonmagnetic spacers: 
near Tc the lattice spins at those interfaces are more 
strongly disordered than the bulk lattice spins, they 
therefore enhance the resistance. The importance 
of this contribution to the enhancement of the total 
resistance depends on the ratio of interface spins to 
bulk spins. This ratio becomes smaller when the 
magnetic layer thickness is larger. 

Far above Tc, most clusters have a small size, so that 
the resistivity is still quite large with respect to the low-T 
phase. However, a few facts account for the decrease of 
the resistivity as T is increased: (i) thermal fluctuations 
are now sufficient to help the itinerant spin overcome en- 
ergy barriers that may occur when it bumps into islands 
of opposite orientation; (ii) the cluster size is now compa- 
rable with the radius Di of the interaction sphere, which 
in turns reduces the height of potential energy barriers. 

We have pitted this interpretation by first creating an 
artificial structure of alternate clusters of opposite spins 
and then injecting itinerant spins into the structure. We 
observed that itinerant spins do advance indeed more 
slowly than in the completely disordered phase (high- 
T paramagnetic phase). This finding is very interesting. 
We believe that it will have other important related phys- 
ical effects yet to be discovered. 

In order to show the existence of clusters at a given 
temperature, we have used the Kopclman algorithm to 
construct clusters^^. We show in Fig. [T2]the distribution 
of cluster size at different temperatures. As can be seen, 
the distribution peak is enlarged with increasing T. 
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FIG. 12: Distribution P{L) of cluster size L at several tem- 
peratures T: (a) above Tc, (b) below Tc. 

We plot in Fig. [T5f a) the cluster size A as a function of 
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FIG. 13: a) Average cluster size versus T b) Average cluster 
size versus (Tc — T) in the In — In scale. 



T. Figure [TSFb) shows the In — In scale of ^ as a function 
of Tc - T. The slope is 0.094 indicating that A does not 
depend significantly on T at T < as expected. 

At this stage, it is worth to make a comparison with a 
recent experiment performed on Ni8iFeig/Al203 multi- 
layer by Brucas and Hanson This system is a magnetic 
permalloy/insulator multilayer which is very similar to 
our model: magnetic layers of thickness t (in the authors' 
notation) are separated by insulator of fixed thickness of 
16 A. Measures of magnetization and resistance have 
been carried out as functions of T for two thicknesses 
t = 16 and 10 A. 

For i = 16 A. X -ray reflectivity, transmission electron 
microscopy and Kerr measurements have shown that the 
magnetic layers are ferromagnetic with Tc ~ 225 K. They 
found that (see Fig. 2a of Ref. [T6h the resistance is very 
small at low T (except for very low T), increases slowly 
with increasing T, makes a jump at Tc and saturated at 
room temperature 300 K. This behavior is very similar to 
what we obtained in the present paper. We note however 
that our model gives a sharper jump than experimen- 
tal data. This is due to the perfect crystalline structure 
(without impurities and defects) of our model which is 
certainly not the case of the experimental system. Be- 
sides, at very low T (< 25 K), due to thermally frozen 
dynamics, experimental measures show an anomaly also 
very similar to MC results at very low T: the decrease 
of A with decreasing T shown in Fig. [5] at T<4 means 
an increase of R with decreasing T. Both experimental 
and theoretical systems show therefore a long-relaxation 
effect due to finite observation time. 
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FIG. 14: Staggered magnetization versus T for several B. 




FIG. 15; Mean free path versus T for several values of B. See 
text for comments 



For t = 10 A, the magnetic layers of the experimental 
system are in fact composed of superparamagnetic do- 
mains. In contrast to the case of t = 16 A, the resistance 
in the case of t = 10 A decreases with increasing T (see 
Fig. 2b of Ref. [TgI ). It is interesting to note that the 
experimental system in this case, which is composed of 
superparamagnetic domains, is equivalent to our model 
in the paramagnetic region above Tc where the existence 
of domains of clusters is shown above. The behavior of 
the resistance observed in our model for T > Tc is in 
excellent agreement with experimental data (see Fig. [TT] 
at T > Tc). The effect of domains on the resistance 
discovered in our present model is thus verified by this 
experiment. 

Finally we show the effect of a magnetic field B applied 
in the z direction. If the inter magnetic film coupling is 
Ji = I, then in the ground state, we need a critical field 
_Bc = 2 to align all spins in the z direction. We show in 
Fig. [14] the lattice staggered magnetization at B=0, 0.5 
and 2. As seen, for B = 2 all lattice spins are aligned in 
the z direction at low T: the staggered magnetization is 
then equal to 1/3 at T = 0. 

An applied field much smaller than Be is expected not 
to modify significantly the itinerant spin current at T <C 
Tc. 

In order to show the effect of the magnetic field 
strength, we define the following quantity termed as 
"magnetoresistance" hereafter 
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FIG. 16: Magnetoresistance versus T for a low field B — 0.5. 
Note the peak at T shghtly larger than Tc. 
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FIG. 17: Magnetoresistance versus T for a large field B = 2. 
Note the peak is larger and higher than that observed in Fig. 

iia 



We now show in Fig. [16] the magnetoresistance in a 
weak field as a function of T. At low T, no significant 
magnetoresistance is expected since itinerant spins are 
parallel to lattice spins. The same is observed at T much 
larger than T^: the lattice spins are parallel to the ap- 
plied field, so itinerant spins will go through the lattice 
without resistance. However, at T slightly larger than Tc 
we observe a large peak of the resistance. This peak has 
the same origin as that observed in Fig. [TTJ namely it is 
due to the existence of the structure of domains in the 
transition region. 

For large fields, the same is observed except that the 
peak is wider and stronger as seen in Fig. [16] for B = 2. 



IV. CONCLUDING REMARKS 

We have studied, by means of MC simulations, the 
transport of itinerant spins interacting with localized lat- 
tice spins in a trilayer system of FCC lattice structure in 
the CIP configuration. Various interactions have been 
taken into account. We found that the spin current is 
strongly dependent on the lattice spin ordering: at low 
T itinerant spins whose direction is parallel (antiparal- 
lel) to the lattice spins yield a strong (weak) current. At 
high temperatures, the lattice spins are disordered, the 
current of itinerant spins is very weak and does not de- 
pend on the input orientation of itinerant spins. As a 
consequence, the resistance is very high at high T. We 
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would like to emphasize here a striking effect found in the 
transition region between low-T ordered phase and high- 
T paramagnetic phase: at T slightly higher than Tc, we 
discovered a peak of the resistance due to the existence 
of domains of lattice spins. Such an existence of domains 
in the critical temperature region is well known from the 
theory of critical phenomena, but no one has expected 
that it would play an important role in the spin transport. 
While writing this paper, we discovered a just-appeared 
experimental paper— which supports our finding on the 
effect of domains in the resistance behavior. 

We have also investigated the effects on the spin cur- 
rent of different parameters which enter in our model: 
nonmagnetic spacer thickness, interaction range, electric 



field, and magnetic field. Physical meaning of each of 
these effects has been discussed. Let us note that so 
far, except Ref. [IB, most magnetoresistance experiments 
have been performed as a function of an applied mag- 
netic field, at a given temperature. While, in our present 
study, we have considered the effect of the lattice order- 
ing on the spin current. We think that in the light of 
the results obtained here, more experiments should be 
performed to investigate the effect of multilayer ordering 
on the spin transport. As a final remark, we note that 
the CPP case is perhaps more difficult to study because 
effects from non magnetic spacers as well as from impu- 
rities and roughness at interfaces will fully set in. Work 
is now in progress to study that case. 
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